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Abstract. "Dual composition", a new method of constructing energy-preserving dis- 
cretizations of conservative PDEs, is introduced. It extends the summation-by-parts 
approach to arbitrary differential operators and conserved quantities. Links to pseu- 
dospectral, Galerkin, antialiasing, and Hamiltonian methods are discussed. 



1. Introduction. For all u,v e C\[-l, 1]), 
\ / vd x w dx — — I wd x v dx + [vw] 1 ^, 



so the operator d x is skew-adjoint on {v G 1,1] : f(±l) = 0} with respect to the 

L? inner product (,). Take n points Xi, a real function v{x), and estimate v'(xi) from 
the values Vi := v (xj). In vector notation, v' = Dv, where D is a differentiation matrix. 
Suppose that the differentiation matrix has the form D = S~ 1 A, in which S induces a 
discrete approximation 

>. 

^ ■ (v, w) s := v T Sw ~ / vw dx = (v, w) 

^ ■ of the inner product. Then 



,T )w, 



(1) (v, Dw) s + (Dv, w) s = v T SS- 1 Aw + v T A T S- T Sw = v T (A + A 1 

which is zero if A is antisymmetric (so that D is skew-adjoint with respect to (,) s ), or 
equals [uiwJLj if x\ = —1, x n = 1, and A + A T is zero except for A nn = —An = |. Eq. ([[]) 
is known as a "summation by parts" formula; it affects the energy flux of methods built 
from D. More generally, preserving structural features such as skew-adjointness leads to 
natural and robust methods. 

Although factorizations D = S~ A are ubiquitous in finite element methods, they have 
been less studied elsewhere. They were introduced for finite difference methods in (see 
0] for more recent developments) and for spectral methods in , in which the connection 
between spectral collocation and Galerkin methods was used to explain the skew-adjoint 
structure of some differentiation matrices. 

Let TC(u) be a continuum conserved quantity, the energy. We consider PDEs 

(2) il = V ^^ 

and corresponding "linear-gradient" spatial discretizations |5], |6], 0, ODEs of the form 

(3) u = L(u)V#(u) 

with appropriate discretizations of u, V, Ti, and 5/5u. For a PDE of the form (0), if 
T>(u) is formally skew-adjoint, then dTi/dt depends only on the total energy flux through 
the boundary; if this flux is zero, 7i is an integral. Analogously, if (|3|) holds, then H = 
UVH) T (L + L T )VH, so that H cannot increase if the symmetric part of L is negative 



definite, and H is an integral if L is antisymmetric. Conversely, all systems with an integral 
can be written in "skew-gradient" form ((|3]) with L antisymmetric) J7|. Hamiltonian 
systems are naturally in the form (|2|) and provide examples. 

This paper summarizes H, which contains proofs and further examples. 



2. Discretizing conservative PDEs. In (Q), we want to allow constant operators such 
as T> = 9" and V = ( _?i o ), and nonconstant ones such as V(u) = ud x + d x u. These differ 
in the class of functions and boundary conditions which make them skew-adjoint, which 
suggests Defn. 1 below. 

Let (, )) be an inner product space. We use two subspaces Jo and #i which can be 
infinite dimensional (in defining a PDE) or finite dimensional (in defining a discretization). 
We write {/_,•} for a basis of $ , {gj} for a basis of 3i, and expand u = Ujfj, collecting 
the coefficients (uj) into a vector u. A cardinal basis is one in which fj(xi) = 5ij, so that 

Uj = u(Xj). 

Definition 1. A linear operator 

V ■ So x $i — > 3, T>(u)v i ^ w, 

is formally skew-adjoint if there is a functional b(u, v, w), depending only on the boundary 
values of u, v, and w and their derivatives up to a finite order, such that 

(v, V{u)w) = — (w, V(u)v) + b(u, v,w) VkG^o, Vw,«)6 #1- 

Ji is called a domain of interior skewness ofV. If b(u,v,w) = Wu G $o, \/v,w G 3i, 
$i is called a domain of skewness ofT>, and we say that T> is skew-adjoint. 

Example 1. Let $ pp (n,r) = {u G C r ([— 1,1]) : u\[ Xi:Xt+1 ] G ^3 n } be the piecewise poly- 
nomials of degree n with r derivatives. For T> = d x , $ pp (n,r), n, r > 0, is a domain of 
interior skewness, i.e., continuity suffices, and {u G $ pp (n,r) : w(±l) = 0} is a domain of 
skewness. 

Example 2. With D(u) = 2(ud x + d x u) + d xxx , we have 

(v, V{u)w) + (w, V{u)v) = [w xx v - w x v x + wv xx + 2uvw], 

so suitable domains of interior skewness are 3o = -5 PP (1,0), #1 = S pp (3,2), i.e., more 
smoothness is required from v and w than from u. A boundary condition which makes 
V{u) skew is {v : v(±l) = 0, ^(1) = ^(-1)}. 

Definition 2. 5o is natural for Ti ifVuE Jo there exists ^ such that 

.. H{u + £v)-H{u) I 6H\ w 

Inn — - K —L = (v,— ) \/veS. 

e^O £ \ dU I 

The naturality of Jo often follows from the vanishing of the boundary terms, if any, 
which appear of the first variation of 7i, together with mild smoothness assumptions. 

We use appropriate spaces Jo and Ji to generate spectral, pseudospectral, and finite 
element discretizations which have discrete energy H := Ti\$ as a conserved quantity. 
The discretization of the differential operator T> is a linear operator T> : Ji — > Jo> and 
the discretization of the variational derivative 4^ is ^ G Ji. Each of D and ^ is a 
weighted residual approximation ||, but each uses spaces of weight functions different 
from its space of trial functions. 



Definition 3. S is the matrix of (, ) |# oX 3i> i- e - <% '■= (fi,9j)- A(u) is the matrix of the 
linear operator A : (v,w) i— > (v,T>(u)w), i.e. Aij(u) := (gi,V(u)gj) . 

Proposition 1. Let Jo be natural for 7i and let S be nonsingular. Then for every u G $0 
there is a unique element ^P- G §1 such that 

Its coordinate representation is S^VH where H(u) := r H.{u i f^). 

Proposition 2. Let S be nonsingular. For every v G Ji ; there exists a unique element 
Vv G Jo satisfying 

(Vv, w) = (Vv, w ) Vi/; G Ji. 
The map v 1— > Vv is linear, with matrix representation D := S~ T A. 

Definition 4. Py^ : Jo ~ > Jo is the dual composition discretization ofV^-. 

Its matrix representation is S'~ T AS' _1 VH. The name "dual composition" comes from 
the dual roles played by #0 an d Ji in defining V and which is necessary so that their 
composition has the required linear-gradient structure. Implementation and accuracy of 
dual composition and Galerkin discretizations are similar. Because they coincide in simple 
cases, such methods are widely used already. 

Proposition 3. // Ji is a domain of skewness, the matrix S^AS' 1 is antisymmetric, 
and the system of ODEs 

(4) ii = S^AS^VH 

has H as an integral. If, in addition, V is constant — i.e., does not depend on u — then the 
system (^) is Hamiltonian. 

The method of dual compositions also yields discretizations of linear differential oper- 
ators V (by taking 7i = | (u,u)), and discretizations of variational derivatives (by taking 
V = 1). It also applies to formally se(/-adjoint V's and to mixed (e.g. advection-diffusion) 
operators, where preserving symmetry gives control of the energy. 

The composition of two weighted residual discretizations is not necessarily itself of 
weighted residual type. The simplest case is when $0 — Ji and we compare the dual 
composition to the Galerkin discretization, a weighted residual discretization of V^- with 
trial functions and weights both in J . They are the same when projecting to Jo, 
applying V, and again projecting to #0, is equivalent to directly projecting Vj^ to Jo. 

For brevity, we assume Jo = Ji for the rest of Section 2. 

Proposition 4. Vj^ is the Galerkin approximation ofV^ if and only ifV(^ — ^) _L 
J . This occurs if (i) £>(5o~) -L -So, or (H) D is exact and applying V and orthogonal 
projection to Jo commute, or (Hi) ^ is exact, i.e., j± G J . 



Fourier spectral methods with T> — d™ satisfy (ii), since then $ has an orthogonal basis 
of eigenf unctions e l ^ x of T>, and differentiating and projecting (dropping the high modes) 
commute. This is illustrated later for the KdV equation. 

The most obvious situation in which ji g #o is when H = | (u,u), since then ^ = «G 
So and = Du, and the discretization of T> is obviously the Galerkin one! When the 
functions fj are nonlocal, D is often called the spectral differentiation matrix. The link 
to standard pseudospectral methods is that some Galerkin methods are pseudospectral. 

Proposition 5. IfV($i) C $ ll then T>v = T>v, i.e., the Galerkin approximation of the 
derivative is exact. If, further, {fj} is a cardinal basis, then D is the standard pseudospec- 
tral differentiation matrix, i.e. D, t j = T>fj(xi). 

We want to emphasize that although A, S, and D depend on the basis, D depends only 
on 3o and 5i, i-e., it is basis and grid independent. In the factorization D = S~ T A, the 
(anti)symmetry of A and S is basis independent, unlike that of D. These points are well 
known in finite elements, less so in pseudospectral methods. 

Example 3 (Fourier differentiation). Let be the trigonometric polynomials of de- 
gree n, which is closed under differentiation (so that Prop. [5]) applies, and is a domain 
of skewness of T> = d x . In any basis, A is antisymmetric. Furthermore, the two popular 
bases, {sin(ja;)™ =1 , cos(ja;)™ =0 }, and the cardinal basis on equally-spaced grid points, are 
both orthogonal, so that S = al and D = S~ 1 A is antisymmetric in both cases. 

Example 4 (Polynomial differentiation). 5i = *}3 n ([— 1,1]) is a domain of interior 
skewness which is closed under T> = d x , so pseudospectral differentiation factors as D = 
S~ 1 A in any basis. For a cardinal basis which includes xq = —1, x n = 1, we have 
(A + A T )ij = —1 for i = j = 0, 1 for i = j = n, and otherwise, making obvious the 
influence of the boundary. For the Chebyshev points Xj = — cos(iir/n), % = 0, . . . , n, A can 
be evaluated first in a basis {Tj} of Chebyshev polynomials: one finds Af^ eh = 2j 2 / (j 2 — i 2 ) 
for i - j odd, and ^ heb - 2{i 2 + j 2 - l)/[((i + j) 2 - l)((i - j) 2 - 1)] for i - j even, 
with other entries 0. Changing to a cardinal basis by Fij = Tj(xi) = cos(ij7r/n), a 
discrete cosine transform, gives A = j?- l j[ cheh j?- T _ p or example, with n = 3 (so that 
(x , xi, x 2 , x 3 ) = (-1, -|, |, 1)), we have 

/ -19 24 -8 3 \ / 4096 -304 496 -1024 \ / -135 184 -72 23 \ 

n _ i 2 -6 -2 6 \ _ c— T 4 _ i / -304 811 -259 496 1 -184 256 -72 \ 

^ — 6 1 -6 2 6 -2 I _ 71 _ 256 1 496 -259 811 -304 J 270 I 72 -256 184 I • 

V -3 8 -24 19 / V -1024 496 -304 4096 / V -23 72 -184 135 / 

S and A may be more amenable to study than D itself. All their eigenvalues are very 
well-behaved; none are spurious. The eigenvalues of A are all imaginary and, as n — » oo, 
uniformly fill [—in, in] (with a single zero eigenvalue corresponding to the Casimir of d x ). 
The eigenvalues of S closely approximate the quadrature weights of the Chebyshev grid. 

For V 7^ d x , V may be quite expensive and no longer pseudospectral. (There is in 
general no S with respect to which the pseudospectral approximation of T>v is skew- 
adjoint.) However, T>v can be computed quickly if fast transforms between cardinal 
and orthonormal bases exist. We evaluate Vv exactly for v G #1 and then project S- 
orthogonally to #1- 

Example 5 (Fast Fourier Galerkin method). Let T>(u) be linear in u, for example, 
T>(u) = ud x + d x u. Let u, v G #1, the trigonometric polynomials of degree n. Then 



T>(u)v is a trigonometric polynomial of degree 2n, the first n modes of which can be 
evaluated exactly using antialiasing and Fourier pseudospectral differentiation. The ap- 
proximation whose error is orthogonal to $1 is just these first n modes, because S = I in 
the spectral basis. That is, the antialiased pseudospectral method is here identical to the 
Galerkin method, and hence skew-adjoint. Antialiasing makes pseudospectral methods 
conservative. This is the case of the linear £>'s of the Euler fluid equations. 



Example 6 (Fast Chebyshev Galerkin method). Let T>{u) be linear in u and let 
u, v G 5i = ^Pn- With respect to the cardinal basis on the Chebyshev grid with n + 1 
points, T>(u)v can be computed in time 0(n\ogn) as follows: (i) Using an FFT, express 
u and v as Chebyshev polynomial series of degree n; (ii) Pad with zeros to get Chebyshev 
polynomial series of formal degree 2n; (iii) Transform back to a Chebyshev grid with 
2n + 1 points; (iv) Compute the pseudospectral approximation of V{u)v on the denser 
grid. Being a polynomial of degree < 2n, the corresponding Chebyshev polynomial series 
is exact; (v) Convert D(u)v to a Legendre polynomial series using a fast transform Pj; 
(vi) Take the first n + 1 terms. This produces V{u)v, because the Legendre polynomials 
are orthogonal, (vii) Convert to a Chebyshev polynomial series with n + 1 terms using a 
fast transform; (viii) Evaluate at the points of the original Chebyshev grid using an FFT. 

3. Examples of the dual composition method. 

Example 7 (The KdV equation), ii + 6uu x + u xxx = with periodic boundary condi- 
tions has features which can be used to illustrate various properties of the dual composition 
method. Consider two of its Hamiltonian forms, 

u = ^i-^p 2?i = d x , Hi = J ( - u 3 + ^u 2 x ) dx, 

and 

u = V 2-^-, v 2 = -{2ud x + 2d x u + d xxx ), H 2 = ^ J u 2 dx. 

In the case 3o — 3l — 5 tng , v := is the orthogonal projection to #o of = 
—3tt 2 — u xx ; this can be computed by multiplying out the Fourier series and dropping all 
but the first n modes, or by antialiasing. Then T>iv = v x , since differentiation is exact 
in $ tng . Since T>\ is constant, the discretization is a Hamiltonian system, and since T>\ is 
exact on constants, it also preserves the Casimir C = J u dx. In this formulation, Prop. |] 
(ii) shows that the dual composition and Galerkin approximations of "Ci^ 1 coincide, for 
differentiation does not map high modes to lower modes, i.e., T>i($ tm&± ) _L J trig . 

In the second Hamiltonian form, H 2 = |u T 5u, = S~ 1 VH2 = u, and the Galerkin 

approximation of is exact, so that Prop. f| (iii) implies that the composition T> 2 ^j L 
also coincides with the Galerkin approximation. V 2 v can evaluated using antialiasing as in 
Example |^. T> 2 is not a Hamiltonian operator, but still generates a skew-gradient system 
with integral H 2 . Thus in this (unusual) case, the Galerkin and antialiased pseudospectral 
methods coincide and have three conserved quantities, Hi, H 2 , and C^ng. 

The situation for finite element methods with $o — -Si — $ pp {n,r) is different. In 
the first form, we need r > 1 to ensure that 3o is natural for Hi, in the second form, 
naturality is no restriction, but we need r > 2 to ensure that is a domain of interior 



skewness. The first dual composition method is still Hamiltonian with integral Hi and 
Casimir C = U{ J fidx, but because D\ does not commute with projection to 3i, it is 
not a standard Galerkin method. In the second form, = u is still exact, so the dual 
composition and Galerkin methods still coincide. However, they are not Hamiltonian. 

Example 8 (An inhomogeneous wave equation). When natural and skew bound- 
ary conditions conflict, it is necessary to take # 7^ $i- Consider q = a(x)p, ft = q xx , 
q x (±l,t) = 0. This is a canonical Hamiltonian system with 




= a(x)p. 



Note that (i) the boundary condition is natural for Tl, and (ii) no boundary conditions 
are required for T> to be skew-adjoint in L 2 . Since j£ is computed with trial functions 
in 3i, we should not include q x (±l) = in #1, for this would be to enforce (— q xx ) x = 0. 
In || we show that a spectrally accurate dual composition method is obtained with 
do = {q e V n+ 2 : q x (±l) = 0} x <p n and ft = ¥n X % n . 

4. Quadrature of Hamiltonians. Computing VH = V7i(ujfj) is not always possible 
in closed form. We would like to approximate 7i itself by quadratures in real space. 
However, even if the discrete H and its gradient are spectrally accurate approximations, 
they cannot always be used to construct spectrally accurate Hamiltonian discretizations. 

In a cardinal basis, let H — J h(u)dx and define the quadrature Hamiltonian H q := 
h(uj)wj = w T /i(u) where Wj = f fjdx are the quadrature weights. Since VH q = 
Wh'(u), ~ W^WHg, Unfortunately, DW~ lJ \7 H q is not a skew-gradient system, while 
DS~ l VH q is skew-gradient, but is not an accurate approximation. 

DW^ 1 VH q can only be a skew-gradient system if DW' 1 is antisymmetric, which occurs 
in three general cases, (i) On a constant grid, W is a multiple of the identity, so if D 
is antisymmetric, DW~ l is too. (ii) On an arbitrary grid with D = (_° 7 o)' DW^ 1 is 
antisymmetric, (iii) On a Legendre grid with # = $1, S = W, and DW~ X = W'^-AW' 1 
is antisymmetric. The required compatibility between D and W remains an intriguing 
and frustrating obstacle to the systematic construction of conservative discretizations of 
strongly nonlinear PDEs. 
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